Humidity-CVD

Data loading

Packages

## data viz for humidity yerawise

pacman::p_load("ggcorrplot","directlabels","pglm", "scales")
suppressPackageStartupMessages(library(tidyverse))

---

Correlation Analysis

# Correlation analysis
combined_data <- combined_data |> 
        rename(maximum_25=high_25,
               maximum_975=high_975,
               minimum_25=low_25,
               minimum_975=low_975,
               average_25=avg_25,
               average_975=avg_975,
               median_25=med_25,
               median_975=med_975,) 
correlation_matrix <- cor(combined_data |> select(average_25, average_975, median_25,median_975,mode_25,mode_975,
                                                  maximum_25,maximum_975,
                                                  minimum_25,minimum_975,variation_25,variation_975), method = "pearson")
#correlation_matrix
# visualize as heatmap


# Compute a matrix of correlation p-values
p.mat <- cor_pmat(correlation_matrix)

ggcorrplot(correlation_matrix, 
           outline.color = "white",
           type = "lower",
           ggtheme = ggplot2::theme_gray,
           colors = c("#6D9EC1", "white", "#E46726"),
           lab = TRUE,
           p.mat = p.mat,
           insig = "blank")

Figure 1: Correlation matrix for humidity variables

This matrix shows pairwise correlations among various humidity metrics (average, median, mode, minimum, maximum, variation) at the 2.5th and 97.5th percentiles. Only statistically significant correlations are displayed, indicating strong associations among several metrics. Based on this analysis, both 2.5 and 97.5 percentiles of the average humidity, and 2.5th percentile of the variation in humidity were selected for further analysis. This is to minimize multicollinearity and retain representative variability in the downstream panel regression.

---

# averages
combined_data |>
        ggplot(aes(x = factor(Year), group=1)) +
        geom_line(aes(y=average_25),linetype="dashed", color="midnightblue")+
        geom_line(aes(y=average_975),linetype="dashed", color="midnightblue")+
        geom_ribbon(aes(ymin=average_25, ymax=average_975),fill="midnightblue", alpha=0.4)+
        facet_wrap(~State)+
        theme_minimal(base_size = 20)+
        labs(y="humidity (%)", x="Year",
             title = "Average humidity",
             subtitle = "2.5th and 97.5th percentiles")+
        theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 12))

# variations
combined_data |>
        ggplot(aes(x = factor(Year), group=1)) +
        geom_line(aes(y=variation_25),linetype="dashed", color="seagreen")+
        geom_line(aes(y=variation_975),linetype="dashed", color="seagreen")+
        geom_ribbon(aes(ymin=variation_25, ymax=variation_975),fill="seagreen", alpha=0.4)+
     #geom_ribbon(aes(ymin=0, ymax=variation_25),fill="seagreen", alpha=0.4)+
        facet_wrap(~State)+
        theme_minimal(base_size = 20)+
        labs(y="humidity (%)", x="Year",
             title = "humidity Variations",
             subtitle = "2.5th percentiles")+
        theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 12))

Figure 2: Humidity summaries for India, 2011-2019

(a) Average humidity levels (2.5th and 97.5th percentiles) across Indian states from 2011 to 2019.

(b) Variations in humidity levels a(2.5th and 97.5th percentiles) across Indian states from 2011 to 2019.

---

model_CVD_Death_avg25 <- pglm(CVD_Death ~ average_25,
                              data = combined_data,
                              family = poisson,
                              model = "within",
                              index = c("State", "Year"))
summary(model_CVD_Death_avg25)
confint(model_CVD_Death_avg25)[, 1]
confint(model_CVD_Death_avg25)[, 2]
irr_avg25_death <- exp(coef(model_CVD_Death_avg25))
irr_avg25_death
irr_lci_avg25_death <- exp(confint(model_CVD_Death_avg25)[, 1])
irr_uci_avg25_death <- exp(confint(model_CVD_Death_avg25)[, 2])
irr_lci_avg25_death
irr_uci_avg25_death

model_CVD_Death_avg975 <- pglm(CVD_Death ~ average_975,
                              data = combined_data,
                              family = poisson,
                              model = "within",
                              index = c("State", "Year"))
summary(model_CVD_Death_avg975)
confint(model_CVD_Death_avg975)[, 1]
confint(model_CVD_Death_avg975)[, 2]
irr_avg975_death <- exp(coef(model_CVD_Death_avg975))
irr_avg975_death
irr_lci_avg975_death <- exp(confint(model_CVD_Death_avg975)[, 1])
irr_uci_avg975_death <- exp(confint(model_CVD_Death_avg975)[, 2])
irr_lci_avg975_death
irr_uci_avg975_death

model_CVD_Death_var25 <- pglm(CVD_Death ~ variation_25,
                               data = combined_data,
                               family = poisson,
                               model = "within",
                               index = c("State", "Year"))
summary(model_CVD_Death_var25)
irr_var25_death <- exp(coef(model_CVD_Death_var25))
irr_var25_death
irr_lci_var25_death <- exp(confint(model_CVD_Death_var25)[, 1])
irr_uci_var25_death <- exp(confint(model_CVD_Death_var25)[, 2])
irr_lci_var25_death
irr_uci_var25_death


### CVD_Incidence

model_CVD_Incidence_avg25 <- pglm(CVD_Incidence ~ average_25,
                                  data = combined_data,
                                  family = poisson,
                                  model = "within",
                                  index = c("State", "Year"))
summary(model_CVD_Incidence_avg25)
irr_avg25_incidence <- exp(coef(model_CVD_Incidence_avg25))
irr_avg25_incidence
irr_lci_avg25_incidence <- exp(confint(model_CVD_Incidence_avg25)[, 1])
irr_uci_avg25_incidence <- exp(confint(model_CVD_Incidence_avg25)[, 2])
irr_lci_avg25_incidence
irr_uci_avg25_incidence



model_CVD_Incidence_avg975 <- pglm(CVD_Incidence ~ average_975,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_Incidence_avg975)
irr_avg975_incidence <- exp(coef(model_CVD_Incidence_avg975))
irr_avg975_incidence
irr_lci_avg975_incidence <- exp(confint(model_CVD_Incidence_avg975)[, 1])
irr_uci_avg975_incidence <- exp(confint(model_CVD_Incidence_avg975)[, 2])
irr_lci_avg975_incidence
irr_uci_avg975_incidence


model_CVD_Incidence_var25 <- pglm(CVD_Incidence ~ variation_25,
                                  data = combined_data,
                                  family = poisson,
                                  model = "within",
                                  index = c("State", "Year"))
summary(model_CVD_Incidence_var25)
irr_var25_incidence <- exp(coef(model_CVD_Incidence_var25))
irr_var25_incidence
irr_lci_var25_incidence <- exp(confint(model_CVD_Incidence_var25)[, 1])
irr_uci_var25_incidence <- exp(confint(model_CVD_Incidence_var25)[, 2])
irr_lci_var25_incidence
irr_uci_var25_incidence


### CVD_Prevalence

model_CVD_Prevalence_avg25 <- pglm(CVD_Prevalence ~ average_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_Prevalence_avg25)
irr_avg25_prevalence <- exp(coef(model_CVD_Prevalence_avg25))
irr_avg25_prevalence
irr_lci_avg25_prevalence <- exp(confint(model_CVD_Prevalence_avg25)[, 1])
irr_uci_avg25_prevalence <- exp(confint(model_CVD_Prevalence_avg25)[, 2])
irr_lci_avg25_prevalence
irr_uci_avg25_prevalence



model_CVD_Prevalence_avg975 <- pglm(CVD_Prevalence ~ average_975,
                                    data = combined_data,
                                    family = poisson,
                                    model = "within",
                                    index = c("State", "Year"))
summary(model_CVD_Prevalence_avg975)
irr_avg975_prevalence <- exp(coef(model_CVD_Prevalence_avg975))
irr_avg975_prevalence
irr_lci_avg975_prevalence <- exp(confint(model_CVD_Prevalence_avg975)[, 1])
irr_uci_avg975_prevalence <- exp(confint(model_CVD_Prevalence_avg975)[, 2])
irr_lci_avg975_prevalence
irr_uci_avg975_prevalence


model_CVD_Prevalence_var25 <- pglm(CVD_Prevalence ~ variation_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_Prevalence_var25)
irr_var25_prevalence <- exp(coef(model_CVD_Prevalence_var25))
irr_var25_prevalence
irr_lci_var25_prevalence <- exp(confint(model_CVD_Prevalence_var25)[, 1])
irr_uci_var25_prevalence <- exp(confint(model_CVD_Prevalence_var25)[, 2])
irr_lci_var25_prevalence
irr_uci_var25_prevalence

IRR - cvd prevalence incidence mortality

Table 1: Incidence Rate Ratios derived from Poisson panel regressions for CVD burden (mortality, prevalence, and incidence)

Humidity	Percentile	Mortality [95% CIs]	Prevalence[95% CIs]	Incidence[95% CIs]
Average	2.5	1.01 [1.004,1.017]	1.016[1.015,1.017]	1.013[1.01,1.016]
Average	97.5	1.002 [1.001,1.003]	1.002[1.002,1.002]	1.002[1.001,1.002]
Variation	2.5	1.014 [0.996,1.032]	1.014[1.01,1.017]	1.011[1.001,1.021]

---

PAF

gbd_data_uncorrected <- read_csv(here::here("data","gbd_data_uncorrected.csv"))
gbd_data_uncorrected <- gbd_data_uncorrected |> 
        filter(Year %in% 2011:2019)

# rename state.name as State
#combined_data_temperature_forpanel 
#<- combined_data_temperature_forpanel |> 
 #dplyr::rename(State = state_name,
 #             Year = year)

# merge temperature data and gbd data

combined_data_uncorrectedgbd <- merge(
  combined_data_relhum_forpanel,
  gbd_data_uncorrected,
  by = c("State", "Year")
)

paf_data_formapviz <- combined_data_uncorrectedgbd |>
dplyr::select(State, Year, CVD_Death, avg_975,
              CVD_Prevalence, CVD_Incidence) |>
filter(Year==2019) |> 
mutate(paf= round((CVD_Death * 0.002)),
       paf_prev= round(CVD_Prevalence * 0.002),
       paf_incid=round(CVD_Incidence * 0.002)) |> 
dplyr::select(-CVD_Death,-avg_975) |> 
        #rename State as state
        dplyr::rename(state = State)

library(sf)
library(tidyverse)

india <- readRDS(here::here("data",
                            "spatial_files", 
                            "India_states.rds"))
india <- india |> 
        rename(state = NAME_1) |> 
        # recode orissa as odisha
        mutate(state = recode(state, "Orissa" = "Odisha")) |> 
# recode uttaranchal as uttarakhand
        mutate(state = recode(state, "Uttaranchal" = "Uttarakhand"))


# combine paf_data_formapviz with india based on state variable
india_paf <- india |> 
  left_join(paf_data_formapviz, by = "state")
india_paf <- india_paf |> 
        mutate(centroid = st_centroid(geometry)) %>% 
  mutate(x = st_coordinates(centroid)[,1], 
         y = st_coordinates(centroid)[,2])


# mortality map
india_paf |> 
  ggplot() +
  geom_sf(aes(fill = paf)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf), size = 3, color = "white") +
  labs(fill = "CVD deaths",
       title = "CVD deaths attributable to high average humidity in Indian states",
       subtitle = "Year 2019") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

# Incidence map
india_paf |> 
  ggplot() +
  geom_sf(aes(fill = paf_incid)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_incid), size = 3, color = "white") +
  labs(fill = "CVD Incidence",
       title = "CVD Incidence attributable to high average humidity in Indian states",
       subtitle = "Year 2019") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())


# Prevalence map
india_paf |> 
  ggplot() +
  geom_sf(aes(fill = paf_prev)) +
  scale_fill_gradient(low = "orange", high = "midnightblue",
                      labels=scales::comma) +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_prev), size = 3, color = "white") +
  labs(fill = "CVD Prevalence",
       title = "CVD Prevalence attributable to high average humidity in Indian states",
       subtitle = "Year 2019") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 3: CVD deaths attributable to high average humidity in Indian states

(a) CVD Mortality

(b) CVD Incidence

(c) CVD Prevalence

---

## daly

model_CVD_DALY_avg25 <- pglm(CVD_DALY ~ average_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_DALY_avg25)
irr_avg25_DALY<- exp(coef(model_CVD_DALY_avg25))
irr_avg25_DALY
irr_lci_avg25_DALY <- exp(confint(model_CVD_DALY_avg25)[, 1])
irr_uci_avg25_DALY <- exp(confint(model_CVD_DALY_avg25)[, 2])
irr_lci_avg25_DALY
irr_uci_avg25_DALY



model_CVD_DALY_avg975 <- pglm(CVD_DALY ~ average_975,
                                    data = combined_data,
                                    family = poisson,
                                    model = "within",
                                    index = c("State", "Year"))
summary(model_CVD_DALY_avg975)
irr_avg975_DALY <- exp(coef(model_CVD_DALY_avg975))
irr_avg975_DALY
irr_lci_avg975_DALY  <- exp(confint(model_CVD_DALY_avg975)[, 1])
irr_uci_avg975_DALY  <- exp(confint(model_CVD_DALY_avg975)[, 2])
irr_lci_avg975_DALY 
irr_uci_avg975_DALY 


model_CVD_DALY_var25 <- pglm(CVD_DALY ~ variation_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_DALY_var25)
irr_var25_DALY <- exp(coef(model_CVD_DALY_var25))
irr_var25_DALY
irr_lci_var25_DALY <- exp(confint(model_CVD_DALY_var25)[, 1])
irr_uci_var25_DALY <- exp(confint(model_CVD_DALY_var25)[, 2])
irr_lci_var25_DALY
irr_uci_var25_DALY


## yld

model_CVD_YLD_avg25 <- pglm(CVD_YLD ~ average_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_YLD_avg25)
irr_avg25_YLD <- exp(coef(model_CVD_YLD_avg25))
irr_avg25_YLD
irr_lci_avg25_YLD <- exp(confint(model_CVD_YLD_avg25)[, 1])
irr_uci_avg25_YLD <- exp(confint(model_CVD_YLD_avg25)[, 2])
irr_lci_avg25_YLD
irr_uci_avg25_YLD



model_CVD_YLD_avg975 <- pglm(CVD_YLD ~ average_975,
                                    data = combined_data,
                                    family = poisson,
                                    model = "within",
                                    index = c("State", "Year"))
summary(model_CVD_YLD_avg975)
irr_avg975_YLD <- exp(coef(model_CVD_YLD_avg975))
irr_avg975_YLD
irr_lci_avg975_YLD  <- exp(confint(model_CVD_YLD_avg975)[, 1])
irr_uci_avg975_YLD  <- exp(confint(model_CVD_YLD_avg975)[, 2])
irr_lci_avg975_YLD 
irr_uci_avg975_YLD 


model_CVD_YLD_var25 <- pglm(CVD_YLD ~ variation_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_YLD_var25)
irr_var25_YLD <- exp(coef(model_CVD_YLD_var25))
irr_var25_YLD
irr_lci_var25_YLD <- exp(confint(model_CVD_YLD_var25)[, 1])
irr_uci_var25_YLD <- exp(confint(model_CVD_YLD_var25)[, 2])
irr_lci_var25_YLD
irr_uci_var25_YLD

## YLL

model_CVD_YLL_avg25 <- pglm(CVD_YLL ~ average_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_YLL_avg25)
irr_avg25_YLL <- exp(coef(model_CVD_YLL_avg25))
irr_avg25_YLL
irr_lci_avg25_YLL <- exp(confint(model_CVD_YLL_avg25)[, 1])
irr_uci_avg25_YLL <- exp(confint(model_CVD_YLL_avg25)[, 2])
irr_lci_avg25_YLL
irr_uci_avg25_YLL



model_CVD_YLL_avg975 <- pglm(CVD_YLL ~ average_975,
                                    data = combined_data,
                                    family = poisson,
                                    model = "within",
                                    index = c("State", "Year"))
summary(model_CVD_YLL_avg975)
irr_avg975_YLL <- exp(coef(model_CVD_YLL_avg975))
irr_avg975_YLL
irr_lci_avg975_YLL  <- exp(confint(model_CVD_YLL_avg975)[, 1])
irr_uci_avg975_YLL  <- exp(confint(model_CVD_YLL_avg975)[, 2])
irr_lci_avg975_YLL 
irr_uci_avg975_YLL 


model_CVD_YLL_var25 <- pglm(CVD_YLL ~ variation_25,
                                   data = combined_data,
                                   family = poisson,
                                   model = "within",
                                   index = c("State", "Year"))
summary(model_CVD_YLL_var25)
irr_var25_YLL <- exp(coef(model_CVD_YLL_var25))
irr_var25_YLL
irr_lci_var25_YLL <- exp(confint(model_CVD_YLL_var25)[, 1])
irr_uci_var25_YLL <- exp(confint(model_CVD_YLL_var25)[, 2])
irr_lci_var25_YLL
irr_uci_var25_YLL

IRR - cvd DALY YLD YLL

Table 2: Incidence Rate Ratios derived from Poisson panel regressions for CVD burden (DALY, YLD, and YLL)

Humidity	Percentile	DALY [95% CIs]	YLD[95% CIs]	YLL[95% CIs]
Average	2.5	1.001 [1,1.003]	1.018[1.012,1.023]	1.001[0.999,1.002]
Average	97.5	1.001 [1.001,1.001]	1.002[1.001,1.003]	1.001[1.001,1.001]
Variation	2.5	1.004 [1,1.007]	1.015[0.999,1.03]	1.003[0.999,1.007]

---

Projection maps for 97.5th percentile of average relative humidity

CVD Deaths

all_forecasts_deaths <- read_csv(here::here("output", "cvd_deaths","combined_forecasts_all_states.csv"))

all_forecasts_2030_deaths <- all_forecasts_deaths |>
  dplyr::filter(Year == 2030) |>
  mutate(paf= round((Forecasted_CVD_Deaths * 0.002),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_deaths)

india <- readRDS(here::here("data",
                            "spatial_files",
                            "India_states.rds"))
india <- india |>
  rename(state = NAME_1) |>
  # recode orissa as odisha
  mutate(state = recode(state, "Orissa" = "Odisha")) |>
  # recode uttaranchal as uttarakhand
  mutate(state = recode(state, "Uttaranchal" = "Uttarakhand"))


# combine paf_data_formapviz with india based on state variable
india_paf_deaths <- india |>
  left_join(all_forecasts_2030_deaths, by = "state")
india_paf_deaths <- india_paf_deaths |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# mortality map
india_paf_deaths |>
  ggplot() +
  geom_sf(aes(fill = paf)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf), size = 3, color = "white") +
  labs(fill = "CVD deaths",
       title = "CVD deaths attributable to high average relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 4

CVD Prevalence

all_forecasts_prev <- read_csv(here::here("output", "cvd_prev","combined_forecasts_all_states.csv"))

all_forecasts_2030_prev <- all_forecasts_prev |>
  dplyr::filter(Year == 2030) |>
  mutate(paf= round((Forecasted_CVD_Prevalences * 0.002),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_prev)

# combine paf_data_formapviz with india based on state variable
india_paf_prev <- india |>
  left_join(all_forecasts_2030_prev, by = "state")
india_paf_prev <- india_paf_prev |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# prevalence map
india_paf_prev |>
  ggplot() +
  geom_sf(aes(fill = paf)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf), size = 3, color = "white") +
  labs(fill = "CVD prevalence",
       title = "CVD prevalence attributable to high average relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 5

CVD Incidence

all_forecasts_incidence <- read_csv(here::here("output", "cvd_inci","combined_forecasts_all_states.csv"))

all_forecasts_2030_incidence <- all_forecasts_incidence |>
  dplyr::filter(Year == 2030) |>
  mutate(paf= round((Forecasted_CVD_Incidences * 0.002),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_incidence)

# combine paf_data_formapviz with india based on state variable
india_paf_inci <- india |>
  left_join(all_forecasts_2030_incidence, by = "state")
india_paf_inci <- india_paf_inci |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# incidence map
india_paf_inci |>
  ggplot() +
  geom_sf(aes(fill = paf)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf), size = 3, color = "white") +
  labs(fill = "CVD incidence",
       title = "CVD incidence attributable to high average relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 6

Projection maps for 2.5th percentile of average relative humidity

CVD Deaths

all_forecasts_deaths <- read_csv(here::here("output", "cvd_deaths","combined_forecasts_all_states.csv"))

all_forecasts_2030_deaths <- all_forecasts_deaths |>
  dplyr::filter(Year == 2030) |>
  mutate(paf_deaths= round((Forecasted_CVD_Deaths * 0.0099),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_deaths)

india <- readRDS(here::here("data",
                            "spatial_files",
                            "India_states.rds"))
india <- india |>
  rename(state = NAME_1) |>
  # recode orissa as odisha
  mutate(state = recode(state, "Orissa" = "Odisha")) |>
  # recode uttaranchal as uttarakhand
  mutate(state = recode(state, "Uttaranchal" = "Uttarakhand"))


# combine paf_data_formapviz with india based on state variable
india_paf_deaths <- india |>
  left_join(all_forecasts_2030_deaths, by = "state")
india_paf_deaths <- india_paf_deaths |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# mortality map
india_paf_deaths |>
  ggplot() +
  geom_sf(aes(fill = paf_deaths)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_deaths), size = 3, color = "white") +
  labs(fill = "CVD deaths",
       title = "CVD deaths attributable to low average relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 7

CVD Prevalence

all_forecasts_prev <- read_csv(here::here("output", "cvd_prev","combined_forecasts_all_states.csv"))

all_forecasts_2030_prev <- all_forecasts_prev |>
  dplyr::filter(Year == 2030) |>
  mutate(paf_prev= round((Forecasted_CVD_Prevalences * 0.0157),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_prev)

# combine paf_data_formapviz with india based on state variable
india_paf_prev <- india |>
  left_join(all_forecasts_2030_prev, by = "state")
india_paf_prev <- india_paf_prev |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# prevalence map
india_paf_prev |>
  ggplot() +
  geom_sf(aes(fill = paf_prev)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_prev), size = 3, color = "white") +
  labs(fill = "CVD prevalence",
       title = "CVD prevalence attributable to low average relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 8

CVD Incidence

all_forecasts_incidence <- read_csv(here::here("output", "cvd_inci","combined_forecasts_all_states.csv"))

all_forecasts_2030_incidence <- all_forecasts_incidence |>
  dplyr::filter(Year == 2030) |>
  mutate(paf_inci= round((Forecasted_CVD_Incidences * 0.0128),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_incidence)

# combine paf_data_formapviz with india based on state variable
india_paf_inci <- india |>
  left_join(all_forecasts_2030_incidence, by = "state")
india_paf_inci <- india_paf_inci |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# incidence map
india_paf_inci |>
  ggplot() +
  geom_sf(aes(fill = paf_inci)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_inci), size = 3, color = "white") +
  labs(fill = "CVD incidence",
       title = "CVD incidence attributable to low average relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 9

Projection maps for 2.5th percentile of variation in relative humidity

CVD Deaths

all_forecasts_deaths <- read_csv(here::here("output", "cvd_deaths","combined_forecasts_all_states.csv"))

all_forecasts_2030_deaths <- all_forecasts_deaths |>
  dplyr::filter(Year == 2030) |>
  mutate(paf_deaths= round((Forecasted_CVD_Deaths * 0.0138),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_deaths)

india <- readRDS(here::here("data",
                            "spatial_files",
                            "India_states.rds"))
india <- india |>
  rename(state = NAME_1) |>
  # recode orissa as odisha
  mutate(state = recode(state, "Orissa" = "Odisha")) |>
  # recode uttaranchal as uttarakhand
  mutate(state = recode(state, "Uttaranchal" = "Uttarakhand"))


# combine paf_data_formapviz with india based on state variable
india_paf_deaths <- india |>
  left_join(all_forecasts_2030_deaths, by = "state")
india_paf_deaths <- india_paf_deaths |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# mortality map
india_paf_deaths |>
  ggplot() +
  geom_sf(aes(fill = paf_deaths)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_deaths), size = 3, color = "white") +
  labs(fill = "CVD deaths",
       title = "CVD deaths attributable to low variation in relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 10

CVD Prevalence

all_forecasts_prev <- read_csv(here::here("output", "cvd_prev","combined_forecasts_all_states.csv"))

all_forecasts_2030_prev <- all_forecasts_prev |>
  dplyr::filter(Year == 2030) |>
  mutate(paf_prev= round((Forecasted_CVD_Prevalences * 0.0138),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_prev)

# combine paf_data_formapviz with india based on state variable
india_paf_prev <- india |>
  left_join(all_forecasts_2030_prev, by = "state")
india_paf_prev <- india_paf_prev |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# prevalence map
india_paf_prev |>
  ggplot() +
  geom_sf(aes(fill = paf_prev)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_prev), size = 3, color = "white") +
  labs(fill = "CVD prevalence",
       title = "CVD prevalence attributable to low variation in relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 11

CVD Incidence

all_forecasts_incidence <- read_csv(here::here("output", "cvd_inci","combined_forecasts_all_states.csv"))

all_forecasts_2030_incidence <- all_forecasts_incidence |>
  dplyr::filter(Year == 2030) |>
  mutate(paf_inci= round((Forecasted_CVD_Incidences * 0.01088),0)) |>
  dplyr::rename(state = State)

rmarkdown::paged_table(all_forecasts_2030_incidence)

# combine paf_data_formapviz with india based on state variable
india_paf_inci <- india |>
  left_join(all_forecasts_2030_incidence, by = "state")
india_paf_inci <- india_paf_inci |>
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[,1],
         y = st_coordinates(centroid)[,2])



# incidence map
india_paf_inci |>
  ggplot() +
  geom_sf(aes(fill = paf_inci)) +
  scale_fill_gradient(low = "orange", high = "midnightblue") +
  # Add paf on the states as text with x and y coordinates
  geom_text(aes(x = x, y = y, label = paf_inci), size = 3, color = "white") +
  labs(fill = "CVD incidence",
       title = "CVD incidence attributable to low variation in relative humidity in Indian states",
       subtitle = "Year 2030") +
  theme_classic() +
  # Make x-axis and y-axis blank
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank())

Figure 12